home *** CD-ROM | disk | FTP | other *** search
/ PsL Monthly 1993 December / PSL Monthly Shareware CD-ROM (December 1993).iso / prgmming / dos / c / snippet.exe / RATKO.C < prev    next >
Encoding:
C/C++ Source or Header  |  1990-06-09  |  8.6 KB  |  192 lines
  1. /*
  2. 224/290 03 Jun 90  22:09:20
  3. From:   Ratko Tomic
  4. To:     All
  5. Subj:   SIEVE 5X8 + PICTURE
  6. Attr:   
  7. ------------------------------------------------
  8. */
  9.  
  10.  /********** PART 1/3: Macros, Globals and count_all() ******************/
  11.  
  12.  /************************************************************************
  13.            SIEVE5x8  (Primes to 1.8 Million in < 60K)
  14.                    Author: Ratko V. Tomic
  15.  *************************************************************************
  16.   This program computes primes using Sieve with bits used for flags and
  17.   by removing multiples of 2,3 and 5 ahead of time. This allows computation
  18.   of primes up to 30*(Number of bytes available). Thus even on a PC with
  19.   520K array (in huge model) you can compute primes up to 15.97 Million,
  20.   good for factoring of numbers up 2.55e14.
  21.  *************************************************************************/
  22.  #include <stdio.h>
  23.  
  24. unsigned _stack = 1024;
  25.  
  26.  typedef unsigned char byte;
  27.  typedef unsigned word;
  28.  
  29.  #define S  60000    /* S*30 = 1,800,000 primes, Small model & 4K stack */
  30.                      /* To increase the limit S, reduce stack size.     */
  31.  
  32.  #define M  (S*8)    /* Maximum bits available                          */
  33.  #define T  (S*30)   /* Upper limit for #'s                             */
  34.  byte a[S+1];        /* Bit map holding S*30 primes, must be all 0's    */
  35.                      /* If iterated, re-initialization is NOT needed.   */
  36.  
  37.  byte x[8]={1,7,11,13,17,19,23,29};   /* Numbers left from 1-30    */
  38.  byte y[8]={12,18,21,27,30,36,47,49}; /* Used to calculate steps   */
  39.  
  40.  #define M6(x)  (((x)+(x)+(x))<<1)              /* Multiply by 6   */
  41.  #define L3(x)  ((unsigned long)(x)<<1<<1<<1)   /* Shift left by 3 */
  42.  #define R3(x)  ((unsigned long)(x)>>1>>1>>1)
  43.  #define R2(x)  ((unsigned long)(x)>>1>>1)
  44.  
  45.  #define cvstep(zz,z,bz) z=((word)R3(zz)), bz=((word)zz&7)
  46.  #define set(j,bj) a[j]|=1<<bj
  47.  #define adv(j,bj,k,bk) j+=k; bj+=bk; if (bj>7) j++, bj-=8
  48.  
  49.  #define bn(n) (n&15?(n&3?(n&1?0:1):(n&4?2:3)):(n&48?(n&16?4:5):(n&64?6:7)))
  50.  
  51.  /*** Count primes after sieve5x8() has finished, called from main()   ***/
  52.  
  53.  long count_all(void)            /* Compute total number of primes found */
  54.  { int c;                        /* which is: 3 + Num_of_0_bits_in_a[]   */
  55.    word w,i;
  56.    long t=3;                       /* Account for primes 2,3,5       */
  57.       for(i=0; i<S/2; i++)         /* Process a[] in word chunks     */
  58.          {                         /* Byte order doesn't matter      */
  59.          w=~*((word*)a+i);         /* There are less 0 bits -> use ~ */
  60.          for (c=0; w; ++c) w&=w-1; /* Count bits set to 1 in w       */
  61.          t+=c;                     /* Accumm. total 0 bits from a[]  */
  62.          }
  63.       return t;
  64.  }
  65.  
  66.  
  67.  /*** PART 2/3: sieve5x8() and utility gen_up2() ****/
  68.  
  69.  sieve5x8()                 /** COMPUTE PRIMES < 30*S **/
  70.  { long P,n,d;
  71.    register word i,bi;
  72.    word j,bm,sav_i,sav_bi;
  73.    static long ss[8];       /* 8 Steps used to cross out multiples of P  */
  74.    static word s[8],b[8];   /* Steps in byte offsets and bit shifts      */
  75.       for (i=0,a[i--]=1;;)
  76.          {
  77.          do i++; while(a[i]==0xFF);       /* Find next 0 bit in a[] */
  78.          if (i>=S) return;                /* End of array reached   */
  79.          bi=~a[i]; sav_i=i;
  80.          bm=-bi&bi;                       /* Bit mask=lowest bit set in bi */
  81.          n=L3(i)+(bi=bn(bi));             /* Get index n for prime found   */
  82.          do {                             /* Process bits in their own loop */
  83.             d=(n-1)&~7L;                  /* == ((n-1)/8)*8;               */
  84.             ss[0]=ss[6]=M6(d)+y[(bi-1)&7];/* == ((n-1)/8)*48+y[(n-1)%8];   */
  85.             if ((d=ss[0]+n)>M) return;    /* Done - no more flags to be set */
  86.             P=(long)i*30+x[bi];;          /* Compute prime P from the index */
  87.             ss[7]=n+n+1;                  /* Compute remaining 7 steps     */
  88.             ss[3]=(R2(n+4)&~1L)+P;        /* == ((n+4)/8)*2+P;             */
  89.             ss[2]=ss[4]=i+((1+P)>>1);     /* == (n/8)+(1+P)/2;             */
  90.             ss[1]=ss[5]=R2(n+2)+P;        /* == (n+2)/4+P;                 */
  91.             for (j=3;j<8;j++)             /* Convert steps to byte step s[] */
  92.               cvstep(ss[j],s[j],b[j]);    /* and bit step b[]              */
  93.             s[0]=s[6];s[1]=s[5];s[2]=s[4];/* Clone remaining 3 byte offsets */
  94.             b[0]=b[6];b[1]=b[5];b[2]=b[4];/* same for bit steps            */
  95.             sav_bi=bi; j=i; i=0;
  96.             do {                          /* Tag mult of P in 8 steps ss[j] */
  97.                adv(j,bi,s[i],b[i]);       /* Advance to next postion in a[] */
  98.                set(j,bi);                 /* Set bit for that position     */
  99.                ++i; i&=7;                 /* Advance to next step (i%8)    */
  100.                }
  101.             while((d+=ss[i])<M);
  102.             i=sav_i; bi=sav_bi;
  103.             do n++, bi++, bm<<=1;
  104.               while(bi<8 && a[i]&bm);     /* Get next 0 bit in this byte */
  105.             }
  106.          while(bi<8);
  107.          }
  108.  }
  109.  
  110.  long gen_up2(long h)  /*** Generate primes<=h, where: 6 < h < T ) ***/
  111.  { word i=0,k=1;       /* This func. is not called for benchmark     */
  112.    long c=3,n,b=0;     /* It assumes sieve5x8() was already called   */
  113.    byte m=2;
  114.       while(1)
  115.          {
  116.          n=b+x[k++];             /* n=Next number to check */
  117.          if (n>h) return(c);
  118.          if (!(a[i]&m)) c++;     /* Prime found (its value is n) */
  119.          if (!(m+=m))            /* Advance to next bit */
  120.            {
  121.            m=1; i++; k=0; b+=30;
  122.            if (i>=S) return(c);
  123.            }
  124.          }
  125.  } 
  126.  
  127.  /** PART 3/3: main(), ITER=5 **/
  128.  
  129.  /****** PC Time calculations ******/
  130.  
  131.  #define pc_time (*(long far*)0x46C)     /* Get IBM tick count */
  132.  time_diff(long t,int *s,int *hs)        /* Compute seconds & 1/100's */
  133.  { long dt;
  134.      if ((dt=pc_time-t)<0) dt+=0x1800B0L; /* Adjust for midnight */
  135.      dt*=1000; *s=(int)(dt/18207);
  136.      *hs=(int)(dt%18207)/182;
  137.      if (*hs>100) *hs=99;
  138.  }
  139.  
  140.  /*** TEST SHELL  ***/
  141.  
  142.  #define ITER 5      /* Iterations for benchmark (on 4.77 PC set to 1) */
  143.  
  144.  #define d3(n) (int)((long)(n)/1000000L)        /* Millions            */
  145.  #define d2(n) (int)(((long)(n)%1000000L)/1000) /* Thousands           */
  146.  #define d1(n) (int)((long)(n)%1000)            /* Ones                */
  147.  long primes;                                   /* Total primes found  */
  148.  
  149.  main()
  150.  { long t0;
  151.    int i,sec,s100;
  152.       printf("Computing primes up to %d,%03d,%03d... ",d3(T),d2(T),d1(T));
  153.       t0=pc_time;                /* Timer ticks used for accuracy */
  154.       for (i=0; i<ITER; i++)
  155.         sieve5x8(),              /* Construct primes < T=1,800,000  */
  156.         primes=count_all();      /* Count number of primes found    */
  157.       time_diff(t0,&sec,&s100);  /* Result should be 135,072 primes */
  158.       printf("Found %d,%03d primes.\n",d2(primes),d1(primes));
  159.       printf("TIME: %d.%02d seconds.\n",sec,s100);
  160.  }
  161. /*
  162.  ---------------------- END OF C CODE -----------------------------------
  163.  
  164.  This is new version for the benchmark I posted recently. The main new
  165.  addon is the new algorithm which eliminates primes 2,3,5 (instead of
  166.  2,3 only). This boosts the packing density from 3 to 3.75 numbers/bit,
  167.  or from 24 to 30 numbers per byte. The new algorithm is in the sieve5x8().
  168.  Some macros were rewritten and some code improvements were done, as result
  169.  it takes roughly same time as sieve3x8() per iteration, even though it's
  170.  more complex and generates more primes (1,800,000 for 60,000 byte array).
  171.  For more accuracy (on a 386/25Mhz), 5 ITERATIONS were used:
  172.  
  173.    Compiler          TIME (sec)   CODE (bytes)    EXE SIZE
  174.   -------------------------------------------------------------- 
  175.    TopSpeed 2.02      27.02          1105           8197  (winner again)
  176.    Watcom 7.0         34.16          1164           6950
  177.    MSC 5.0            35.09          1266           8827
  178.    Turbo C 2.0        35.75          1246           8056  (last again)
  179.   --------------------------------------------------------------
  180.  
  181.  THE BIG PICTURE - If you have VGA or MCGA video adapter try this: Take 
  182.  the content of a[S] word at a time (2 bytes) and use it as an offset
  183.  into video RAM, (video mode 0x13: 320x200, 256 colors, flat 2-D address).
  184.  At that offset increment a byte to new (brighter) color. Do that for all
  185.  30,000 words in a[].  Quite surprising ?!  It also works with sieve3x8()
  186.  and other video cards. Also, the byte ordering in a word doesn't matter.
  187.  
  188. --- TBBS v2.1/NM
  189.  * Origin: MSI S/W BBS-Fram. MA-(508) 875-8009-CodeRunneR-HandiWARE (322/327)
  190. */
  191.  
  192.